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Abstract 

We propose a novel approach for the determination of the nature of ultra-high energy cosmic rays by 
exploiting the geomagnetic deviation of muons in nearly horizontal showers. The distribution of the muons 
at ground level is well described by a simple parametrization providing a few shape parameters tightly 
correlated to -^^ax? depth of maximal muon production, which is a mass indicator tightly correlated to 
the usual parameter the depth of maximal development of the shower. We show that some constraints 

can be set on the predictions of hadronic models, especially by combining the geomagnetic distortion with 
standard measurement of the longitudinal profile. We discuss the precision needed to obtain significant 
results, and we propose a schematic layout of a detector. 


1. Introduction 

The nature of the UHECR (cosmic rays of ultra-high energy, of the order of 10^^ eV or more) is one of the 
most challenging and important open questions in astrophysics, as it is crucial for the understanding of their 
origin and of the mechanism of their production. The UHECR are observed through the detection of the 
cascade of particles they induce in the atmosphere: the first interactions occur at energies which cannot be 
reached in colliders, so they are described by models extrapolated above the domain where they can be fitted 
to the data. As a result the mass composition of UHECR is still ambiguous and is motivating considerable 
efforts to improve the observation and the combination of various mass indieators and to reduce the model 
dependent uncertainties. Several experimental results have been obtained so far based on the observation 
of the shower development (e.g., depth of the shower maximum, depth of muon production), the particle 
content and the ground (see for example some recent publications [1013 H 0 El [7]). In this paper we 
propose a novel method for the determination of the mass composition and the hadronic interaction models 
by exploring the features of horizontal air showers, impacting the atmosphere with a zenith angle larger 
than about 60°. 

One key point to identify the primary particle is evaluate the muonic component of the atmospheric 
shower. It is negligible within the core, and has to be evaluated in ground detectors at remote distance 
from the shower axis. At moderate zenith angle, the incident flux on the ground is essentially a mixture 
of photons, electrons, positrons (the eleetromagnetie component) and muons, and it is difficult to separate 
unambiguously the muonic components. At large zenith angle, the slant depth of the atmosphere is large, 
so the electromagnetic component is extinguished at ground level. The aim of this paper is to show how the 
observation of the geomagnetic deviation of muons in nearly horizontal showers may help to reduce these 
uncertainties and provide useful cross-checks between different measurements of mass indicators. 

The paper is organized as follows: In Sect. we describe the atmospheric shower and we explain the 
principle of the method; the approach used to simulate the muonic flux at ground is explained in Sect. 
and we introduce in Sect. a convenient parameterization of the generated maps. In Sect.[^ we show that 
the coefficients of this parametrization are tightly correlated to the nature of the primary particle and we 
examine their dependence on the models of hadronic interactions; in Sect. we evaluate the precision 
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needed on ground measurements to obtain predictive results and, in Sect.[^ we propose possible layouts for 
a detector. 

2. Principle of the method 

Extensive atmospheric showers produced by protons or nuclei include a hadronic component (mainly 
charged mesons) and an electromagnetic component (photons, electrons and positrons) induced by the 
decay of neutral mesons (mainly pious) into photons. These components reach their maximal development 
after crossing between 500 and 1000 g/cm^, depending on the nature and the energy of the primary particle, 
and then decrease progressively. The depth of maximal size Xmax (where the number of charged particles, 
mainly electrons and positrons, is maximal) is known and exploited for a long time, and may be considered 
as a standard; it is essentially sensitive to the electromagnetic component. 

For moderately inclined showers, both components reach ground level, and their relative importance 
(especially the muonic content) is in principle an indicator of the nature of the primary, but the techniques 
needed to evaluate separately the components are not straightforward. At large zenith angles 0 (typically 
6> > 60 deg) the electronic cascade is extinguished and the hadronic one has been transformed through meson 
decays into a flux of muons which do not interact strongly. These muons lose their energy mainly through 
ionization, and many of them decay in flight, but a significant fraction reaches the ground. Their initial 
azimuthal distribution is uniform around the shower axis; but this symmetry is broken by the deviation 
in the magnetic field of the Earth, which results in a distortion of the density at ground level. A detailed 
discussion of this effect may be found in [8]. 

The distribution of the muons in altitude and energy is illustrated in Fig. for showers with a zenith 
angle around 75 deg, the energy of most muons reaching the ground is of the order of a few GeV to a few 10 
GeV, and their path is a few 10 km, so the magnetic deviation is of the order of a few 100 m: the distortion 
is sizeable and may be measured by a surface detector with a sufficient granularity. 

The deviation is proportional to the square of the path of the muons down to the ground and to the 
inverse of their energy, which increases in average with the path because of the losses through decay. The 
net effect is an increase of the distortion with and we can expect the distortion to depend also on the 
longitudinal evolution of the cascade, and especially on the depth of maximal production of muons, 

which is tightly correlated to the standard parameter X^ax foi* a given value of the zenith angle 6>, as it is 
shown in Fig. 

The dependence on 0 is due to a density effect: for higher Xmax is reached at higher altitude, where 
the interaction length is larger, so the mesons decay earlier in terms of atmospheric depth. Both Xmax and 
X^ax indicators of the nature of the primary particle when using a reliable model of the interactions at 
very high energy, especially the hadronic ones. 

Our aim is to express this dependence through a simple paramerization of the muon density and provide 
tools to extract from the observation of inclined showers an estimator of the mass composition using a given 
model of hadronic interactions at ultra-high energy; we want also to obtain some discrimination between 
different models through their prediction of the correlation between X^ax and This approach has the 

advantage of using the pure muonic component of the shower. Moreover in inclined showers the hadronic 
cascade is observed at a higher energy level than in nearly vertical showers, because the mesons decay earlier. 
So the information is complementary to studies at low zenith angles. 

Taking average values of the parameters, we obtain a simple analytic expression of the density that can 
be used as an alternative to previous propositions ( 0 , 191 ), which did not account for the variation of the 
longitudinal profile. 
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Figure 1: Distribution in altitude and energy of the muons produced in inclined showers; in grey: muons reaching the ground. 
Left: proton shower of 1 EeV, 6 = 72 deg; right: 80 deg. 
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Figure 2: Correlation between the depth of maximal size Xmax and the depth o maximal muon production X^ax- Left: Showers 
at ^ = 72 deg, from protons (open symbols) or iron nuclei (solid symbols) using different models for hadronic interactions (see 
below Sect. |5.2| Right: dependence on zenith angle for the QGSJET II model. 


The geometrical frames and coordinates used hereafter are described in Fig. 
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Figure 3: Geometry of an inclined shower: the front plane is perpendicular to the shower axis, which makes an angle 6 with 
the vertical direction. A point in the front plane is defined by the polar coordinates r, 'ijj. 


3. Simulation of the muonic component 
3.1. General procedure 

In this paper we estimate the muon density through a procedure in two stages, to reduce the computing 
load: 

• Extract from a shower simulation package the list of the muons at their production point. We simulate 
proton and iron showers at different primary energies between 10^^ and 10^^ eV, and zenith angles 
between 64 and 80 deg. 

• From each sample, propagate the muons to obtain the density on the ground or in a front plane 
(orthogonal to the shower axis, see Fig. [^, with different values of the transverse magnetic field, 
which is the only relevant quantity because the muons that reach the ground are nearly parallel to the 
axis. 

We do not try to describe the core of the shower (distance less than about 100 m), nor the behaviour at 
large distances (^ 1 km), where the density is well below one muon per square meter, so that detectors of 
a reasonable size have little chance to make a precise measurement. We suppose that the detectors cannot 
distinguish the charge of the muons, so we consider only the total flux as measurable. 

To transpose the geometrical observations into a description in terms of the depth X, the vertical profile 
of the atmosphere as a function of the altitude X{h) needs to be known at any time, not only to define 
the stage of evolution of the hadronic cascade, but also to describe properly the rate of decay of mesons, 
which depends on the density dX/ dh . Fortunately, this profile has little diurnal and seasonal variations at 
altitudes between 10 and 30 km above sea level uni, where most of the muons are produced in the showers 
we are interested to. So in the simulation of the showers we use the standard Linsley atmosphere m- For 
the propagation of the muons we use either the same model, or a pure exponential profile to obtain analytical 
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evaluations. 


Through their decays and their radiative interactions initiating electromagnetic cascades of low energy, 
the muons generates continuously the so-called “electromagnetic halo” (see for example ED- These cascades 
develop over a short distance compared to the muon path to the ground, so the electromagnetic flux follows 
closely the shape of the muonic one, and undergoes the same magnetic distortion. The impact on surface 
measurements (essentially a constant factor in most cases) depend on the nature of the detectors, and is 
beyond the scope of this paper. So we do not simulate the electromagnetic halo. 

The energy loss £ = —dE/dX is supposed to be constant (2 MeV.g“^cm^). This is a good approximation 
in the medium range of energy (few 100 MeV to 10 GeV). Very energetic muons lose more per unit of depth, 
but the total loss is anyway a small fraction of their initial energy; moreover they are concentrated in the 
core, and are weakly deviated. Muons of low energy stop rapidly and contribute little to the density at 
ground. 


3.2. Simulation tools 

To simulate the atmospheric shower we use the package CORSIKA [13]. We set options for the thinning 
procedure such that the statistical weight of the hadrons remains low (less than 10 within the pure hadronic 
cascade and less than 1000 for hadronic subproducts of the electromagnetic cascade through a photo-nuclear 
interaction). 

CORSIKA offers an option to output the list of the muons generated by the hadronic cascade, with their 
production point, their momentum at this point and their statistical weight w inherited from the parent 
particle. The production region of the muons reaching the ground is mainly concentrated within a few ten 
meters around the shower (see Fig. [^, and this distance is small compared to both the magnetic deviation 
and the dispersion due to the multiple Coulomb scattering. So we keep only the longitudinal position of the 
origin of the muon. 
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Figure 4: Distribution in distance from axis vs energy at production point, for muons produced in proton showers of 1 EeV, 
6 = 72 deg, reaching the ground. 


The original direction of the muon is slightly affected by the magnetic deviation of its parent (most of 
the time a pion or a kaon), which has the same charge and a similar energy. For a given momentum p the 
mean flight is I = rpjm; the ratio r/m is much smaller for a kaon than for a pion, so we will discuss only 
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the pion case. The angular deviation over a length I is 6 = cBJlp where Bf is the transverse field (in Tesla) 
and p the momentum in eV/c, that is 6 = cBit jm\ with Bt < 6.10“^T for the field at the surface of the 
Earth, we see that S < 1 mrad anywhere. This value is small compared to the deviation of the produced 
muon along its path to the ground, and we have checked that changing systematically the initial directions 
of the muons by ±1 mrad (depending on the charge) does not modify significantly our results. So we assume 
a perfect symmetry of the directions around the axis and to obtain smooth density maps we replace every 
muon by N clones obtained by successive rotations of 27t/N around the shower axis, each one with a weight 
w/N. 

This sample is the input for two different procedures aimed to obtain the map of the muon fiux on the 
ground: 

• a stepwise extrapolation of the muons through the atmosphere down to the ground, accounting for 
energy loss, multiple scattering, decay probability and magnetic deviation in each step. For a muon 
of weight re, N is chosen as the integer just above re/0.2, so the clones have a weight less than 0.2. 

• an analytic computation of a density of probability in the front plane for each input muon, in the 
approximation of an exponential atmosphere profile, and a weighted summation of these densities. In 
this option N = 24. 

These procedures are detailed in the following subsections. 

3.3. Stepwise approximation 

Each clone is propagated by steps (initially 1 km) until the muon stops, decays or reaches the ground. 
The energy loss is computed according to the local density of air, which is defined as a function of the 
altitude. The step length is divided by 2 if the muon loses more than the half of its energy. In each step 
the position and the direction undergo a deterministic deviation from the magnetic field, and a random one 
to account for the multiple scattering [16]. When a step goes below the ground level, a linear interpolation 
gives the impact at the ground. It was checked that reducing the step length did not modify significantly 
the final distribution. 

From the set of surviving muons, we want to define a ground density. A first option is just to count 
the number of muons per unit of ground surface. Doing so we observe a forward-backward asymmetry due 
to the divergence of the muons: the muons hitting the ground in the upstream region have a large angle 
of incidence. We can also define an intrinsie density as the number of muons hitting a unit surface in a 
plane orthogonal to the shower axis. In that case there is still an asymmetry because the upstream and 
downstream impacts are not at the same longitudinal position in the shower. In practice what matters is 
the number of muons entering a detector; if the detector is a planar vertical surface, the intrinsic density is 
more relevant to estimate the response. We will show in Sect. |4.3| that the forward-backward asymmetry 
does not prevent us from exploiting the magnetic distortion in a reliable way. 

3.4- Analytie ealeulations for an exponential atmosphere 

In the local ground frame {z axis upwards on the vertical direction, with the ground at z = 0), the 
atmospheric depth (quantity of matter above altitude z) is described by : 

X{z)=XaeM-^/L) ( 1 ) 

where L is the attenuation length (typically 8 km). We define the shower frame with a coordinate z' along 
the shower axis (origin at ground level), y' perpendicular to the shower axis in the horizontal plane (see Fig. 
[^. An inclined shower with a zenith angle 6 is considered as equivalent to a vertical one in an exponential 
atmosphere with a scaled attenuation length: 

cosO \ L/cosOJ \ Lsi J 

6 




where Xsi is the slant depth at ground level and Lsi the attenuation length along the axis. Actually there is a 
transverse gradient of density in the shower frame, so this expression is valid only around the axis, at a short 
distance compared to L: we show in Sect. |4.3| that this asymmetry, combined with the fact that the ground 
is replaced by a plane orthogonal to the shower axis, has little consequence on the evaluation of the distortion. 


Let us consider a muon (mass m, lifetime r) injected at that is at a slant depth Xi = Xgi exp(—z'/Lg/)), 
with a kinetic energy Ei^ nearly parallel to the shower axis. We define A = cr jm^ • + eXi (energy 

extrapolated backwards to infinite altitude). 

In Appendix we obtain the following results: 


The muon can reach the ground if Egr = E^o — sXsi > 0, with a probability 


P — 

JTgr — 


EgrXj ^ 

EiX.,^ 


(this expression was already derived in [14] within the same approximations) 


The magnetic angular deviation (perpendicular to the transverse field ^) is: 


^ ~ T- a ( in 1 

Eoo^o^O \ EgrJ 

(with P = ecBt or = cBt if the energies are expressed in eV) 
and the transverse displacement from a straight line: 


(5 = 




E^ 


Fi{a,z'jLsi) with a = , Fi{a,C) = [ — 

-t^oo Jo t 


udu 


(aexp(— 1 ^) 


(3) 


(4) 


(5) 


• The variance of the angular deviation due to the multiple scattering is: 

2 _ El, (^ n 

eXrad \Egr Ei) 

and the variance of the displacement: 

<^pos = -^^F2{a,zJL,i) 

V -^oo / -^rad 

with 

du 


(6) 

(7) 

(8) 


With these expressions we can attribute a weighted density in the {x' ^y') plane to each muon emitted 
with a weight w at z'- along the unit vector {u'^^u'y^u'^); for example, if Bt is along y' axis, so the deviation 
along x': 


f{x',y') =wP, 




exp 


{x'-u',z',±6)^ + {y'-u'^zlf 


where ± represents the sign of the muon charge. 


(9) 
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4. Parametrization of ground densities 
4-1. General features 

As expected, the magnetic distortion increases with the zenith angle and with the magnitude of the 
transverse field Bt. This is illustrated in Fig. which shows the contour levels of the densities of muons 
(total or charge separated). The dependence of Bt is not trivial; we can distinguish two extreme cases: 

• weak distortion (low 0 and Bt)'. the total density is the sum of the positive and the negative com¬ 
ponents: at a given position, each one differs from the undistorted one by a small amount, which is 
proportional to Bt^ with opposite signs, so the global effect cancels at first order, and the variation is 
of the second order in Bt. 

• strong distortion (large 0 and Bt): almost everywhere, the flux is fully dominated by muons of one 
sign, so the dependence on Bt is stronger. 



Figure 5: Contour levels of the muon density in the transverse plane, for a proton shower of 10 EeV at 3 zenith angles (from 
left to right: 64, 72 and 80 deg), with a transverse field of 10 or 60 /rT along y axis. In red: in h\ue:y~, in black (dashed): 

total. The lines correspond to equidistant levels in log scale (2 per decade), starting from 10“^ muons/m^. 


4 . 2 . Funetional parametrization 

For each value of 0 and Bt we want to find an analytical expression of the density in the front plane at 
ground level, as a function of the distance r to axis and the azimuthal angle -0, with 7 /; = 0 in the direction 
perpendicular to Figure and show that for a wide domain in the useful ranges of 0 and Bt the 
dependence of the logarithm of density is approximately linear in yT and sinusoidal in with a relative 
amplitude varying smoothly with r. So, introducing a reference distance r^-ef (typically the average distance 
where the density may be measured) and the variable p = ^Jr/ r^ef — 1, the density is well fitted by: 

/(r, 7 /)) = exp (A(p) + ol{p) cos(2('0 — f^s))) direction of the deviation) (10) 
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and the dependence on p is well decribed by a parabolic parametrization: 

A(/:)) = Aq + Aip + \2p^ 1 <^(p) = olq + OL\p + (y,2p^ 


( 11 ) 


Ao is a size parameter, roughly proportional to the primary energy, while the shape parameters Ai, o^o, 
ai carry most of the information on the shape of the ground density that can be extracted from the 
measurements; a 2 and A 2 are relatively small and need a wide measurement range in r to provide a useful 
information. The values of Ai, q^o, oii are displayed in Fig. [^as functions of 0 and Bt for an average over 
10 proton showers of 1 EeV. 





p 


p 


Figure 6: Parametrization of the muon density in (r, b) coordinates, for a proton shower at 10 EeV, ^ = 64 deg, Bt = 30 p,T. 
Top left: density as a function of ^/r and to p right : azimuthal dependence at different distances; bottom left: parameter 

A (logarithm of the density) as a function oi p = y/r/r^ef — 1; bottom right: parameter a (azimuthal modulation) as a function 
of p. 
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Figure 7: Parametrization of the muon density in coordinates, for a proton shower at 10 EeV, ^ = 80 deg, Bt = 60 /xT 

(same organisation as Fig. |^. 



1 


-3 


O 80 deg ^ 

0.9 

L 

-3.2 


V 76 deg 


: o 


: o 

A 72 deg 

0.8 

- o 

-3.4 

1 o 

■ 68 deg q 




o 

@ 64 deg 

0.7 

7 T 

-3.6 

T 



: o 


T 


0.6 

1 

-3.8 

1 

O 


T 


O T 

T 






0.5 

- O A 

-4 

A 

T 



T 


A 

T 

O 

0.4 

A ■ 

-4.2 

: ° . 

T A 


O 

◄ 

► 

■ 

• 


■ 


0.3 

— 

-4.4 


■ 

◄ 

► 

O 


■ • 

A 


^ • 

A ■ 

• 

0.2 

" , • 

-4.6 

r ^ ■ • 

T A ■ • 


- o ^ ^ 


■ • 

o , ■ • 

0.1 

A # 

■ 

-4.8 

• 

• 

- , I , » , V , , , , , 

0 

7 i , * . 

-5 

. . , . I . T , , 1 T , , 


Bt (a^T) 


Bt (mT) 


Bt (mT) 


Figure 8: Dependence of the parameters of the muon density on the zenith angle 6 and the transverse magnetic field Bt. 
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J^.S. Comparison of the analytical density to the stepwise one 

By construction the analytic density has no forward-backward asymmetry in the shower frame, while the 
stepwise evaluation includes two sources of asymmetry: the non uniform density of air and the fact that the 
ground is not perpendicular to the shower axis. Moreover the distortion depends not only on \Bt\, but also 
on its direction in the front plane. To account for the asymmetry in the stepwise density, we introduce 
in the parametrization a term in cost/^, that is: 

/{r, Ip) = exp (A(/9) + a{p) cos(2(t/) - iPb)) + P(p) cos{^p)) (12) 

In each slice in r we fit A, a and /3 with this function of 'ip. An example is given in Fig. which shows that 
the asymmetry is almost suppressed when going from the ground densities to the intrinsic ones, as defined 
in Sect. |3.3[ that is, the ground asymmetry is dominated by the divergence of the muons from the shower 
axis. For both options we can extract a parabolic parametrization of A and a as functions of p, that is 
coefficients Aq, Ai, A 2 and ao, <ai, 0 ^ 2 , which are quite similar to the ones obtained with the analytic method. 



Figure 9: Density of muons at r = 1000 m for a shower of 1 EeV, at 0 = 72 deg, with Bt = 60 /rT, obtained from the stepwise 
approximation, for different orientations of the transverse magnetic field. Solid: density on ground; open: intrinsic density 
(projected onto the front plane) 


5. Correlation of the shape parameters with the nature of the primary 


5.1. Dependence on shower evolution 

Here we choose QGSJET 11.04 [15] as a reference model for the hadronic interactions. With 0 = 72 deg, 
Bt = 30 juT and Eprim = 0-1, 1 and 10 EeV (10 proton and 10 iron showers at each energy), Eig. 10 shows 
a tight correlation between the relevant shape parameters (oq, oi, Ai) and the depth of maximum 

muon production. It is important to note that proton and iron showers at different energies are on the 
same line, that is, the shape parameters provide an indirect measurement of that can be used for the 

identification of the primary. The same behaviour is observed for the other values of the zenith angle and 
the transverse field. 

Of course the effective identification power (for individual events, or statistically for a sample of events) 
depends on the precision that may be achieved on the measurement of the zenith angle and the shape 
parameters. This will be discussed in Sect. 
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Figure 10: Dependence of the shape parameters on for proton and iron showers at different energies, 6 = 72 deg, for 

the QGSJET 11-04 model. The size of the symbols is related to the energy: small for 0.1 EeV, medium for 1 EeV, large for 10 
EeV. 


In principle the density computed with the stepwise extrapolation is more realistic than the analytic one. 


Fig. 11 shows that the shape parameters of the density at ground level have the same dependence on 


as found with the analytic method. The difference is essentially a global shift, depending on the direction 
of the field in the front plane. So we did not apply the stepwise method to the whole set of showers, and we 
draw conclusions from the results of the analytic method. 


12 







(g/cm^) (g/cm") (g/cm") 



x'^-.ax (g/cm") (g/ctn") x^,,, (g/cm^) 


Figure 11: Depencence of the shape parameters ao, ai, Ai on X^ax, obtained by the stepwise extrapolation for different 
orientations of the transverse field (open symbols), and by the analytic method (solid poins), applied to a proton shower of 1 
EeV, at 72 deg. Top: Bt = 30 p,T; bottom: Bt = 60 fiT. 


5.2. Dependence on hadronic model 

We have chosen 3 available alternative options: EPOS-LHC [17], QGSJETl [T8| and SIBYLL [T9|, the 
last two being generally considered as obsolete. Eor each one, we have simulated 10 proton and iron showers 
at 1 EeV, 72 deg. Eig. shows that for each model the shape parameters have a similar correlation to 
t)ut the lines differ significantly for at least one of the parameters. In some cases, the distribution 
of one parameter may in principle discriminate two models, whatever the composition. Eor example, the 
distribution of ai separates well QGSJETl and QGSJET 11.04 on the one side from EPOS and SIBYLL on 
the other side. Discriminating variables may be defined as combinations of ao, (ai, Ai. Again the effective 
discrimination power depends on the precision of the measurement of these parameters (see Sect. |^. 
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Figure 12: Dependence of the shape parameters on X^ax at F' = 1 EeV, 6 = 72 deg, Bt = 30 /rT, for different hadronic models. 
Blue symbols are for iron, red ones for proton. 


We have also tried to modify by hand some characteristics of the muons. First, by multiplying all 
transverse momenta by 1.1, because the distribution in pt may affect the interpretation of the muonic flux 
at a fixed distance from core. We have also simulated a global contraction of the longitudinal profile of 
muon production, without modifying by replacing for each muon the depth of production Xi by 

^max + 0-9 * {^i — ^max)’ dividing its energy by the ratio of the air density at the new position to the 
original one. The results are shown in Fig. the contraction of the profile has no significant effect, while 
the scaling of pt results in an important shift of the parameters. This means that the magnetic distortion 
is very sensitive to the distribution of the transverse momentum, which is governed by the latest hadronic 
interactions. 







X^n,ox 


X".ox 


Figure 13: Dependence of the shape parameters on X^ax at E = 1 EeV, 6 = 72 deg, Bt = 30 yuT, with QGSJET 11.04, on 
artificial modifications. Solid: no modification; open circles: transverse momentum of muons multiplied by 1.1; open triangles: 
longitudinal profile contracted by 0.9 (see text). 


6. Precision on the measurements 

We want here to obtain an evaluation of the measurement errors using an ideal detector which is a regular 
array of identical elements acting as pure muon counters, that is, giving the number of muons crossing their 
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area A (projected onto the front plane of the shower axis) with a poissonian distribution. Actually the 
detector may see the electrons of the electromagnetic halo, and also some photons if they interact with the 
material: from this point of view, this approximation is conservative. 

For the following evaluations we define a standard detector array: a square array of spacing £ = 500 m on 
a horizontal ground, with A = 10 m^. Actually the precision depends essentially on the density of detectors 
in projection onto the front plane, that is 1 /{£‘^ cos 0 ) for the standard array. 


6.1. Precision on zenith angle and transverse field 

As can be seen in Fig. the shape parameters depend strongly on the zenith angle 6 : if we want to 
make an event-by-event study, we have to be sensitive to a variation of of the order of 30 g/cm^, so we 
need to measure 0 with a sufficient precision: for example, around 72 deg, with Bf = 30 /iT, the variation 
of ao is of the order of 0.04 per degree, which correspond to a variation of 60 g/cm^, so a precision better 
than 0.5 degree fulfills the requirement. At higher values of 0 and the dependence in 0 is much stronger; 
this is partially compensated by the fact that the dependence of ao on is also increased, but not fully. 

Fortunately the shower front becomes thinner and flatter with increasing so the error on the direction 
decreases. 

The precision on 0 may be evaluated for the standard array through a conservative approach: 

• making a precise evaluation of X^^^^ for a single event requires to have measurements on a wide area, 
at least up to 1000 m from shower axis. 

• we assume that the precision on the front time measured in a detector is the dispersion of the arrival 
times of the muons, instead of accounting for the time of the first one, which has less fluctuation in 
case of multiples hits. 

• we ignore the information from detectors at more than 1000 m from the core and we attribute to the 
measurements the dispersion at found at 1000 m by the stepwise propagation of muons. This is an 
overestimation because the dispersion increases with the distance. We find around 100 ns for 0 = 64 
deg, 50 ns for 0 = 72 deg and 30 ns for ^ = 80 deg. We can suppose that the timing precision of the 
detector is better than 10 ns. 


With these values, we can evaluate the precision on 0 using the projection onto the front plane. For example, 
assuming the core to be at the center of one square of the array, with x' axis along the forward-backward 
direction, we keep ^ l/(-^/1000)^ cosO detectors with transverse coordinates x[ = 4z{£/2) cos6>, ±(3-^/2) cos6>, 
etc, and we obtain ag through a summation over them: 


1 _ ^ / x' y _ 2 

al ^ yccTt / ^ 

^ | cc '|<1000 ^ ^ 


3 V cat cos 0 




(13) 


(this error scales as l/£). With the above values for at, we find for example es ao 0.5 deg at 6> = 64 deg, 
0.2 deg dX 0 = 12 deg and 0.1 deg at 6> = 80. This is sufficient for our purpose. 

If a sample of n events is used statistically, what matters is ao /so the requirement of the footprint 
on the ground may be relaxed. In any case systematic errors, for example the bias due to the ground 
asymmetry, should kept be under control. 


6.2. Precision on the shape parameters 

We suppose that the flux seen by the detectors may be computed in the front plane {r^fi) with the 
approximation of Sect. so that the parametrization found there may be used for the average number of 
muons in a detector. 
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Due to the spacing of the detectors and the rapid decrease at large distance, the range in r providing 
effective information is limited; If we choose for r^ei the medium value of ^/r^ we can omit the terms in 0^2 
and A 2 and write the expected number of muons in a detector at position (r, 2 p) as: 


N{r, = exp (Ao + Xip + ( 0^0 + c^ip) cos(2('0 - ^^b))) 


with 



(14) 


The number of muons in a detector at position follows a Poisson law of mean value N{ri^^pi] p) 

depending on the parameters p = (xc, ^ 0 ? <^05 <^ 1)5 where Xc^Vc is the core position. For given values 

of 0 and Bf we can also apply this formalism to a “compound” parametrization where Ai, o^o and ai 
are expressed as linear functions of according the dependence observed in Sect. |5.1| (Fig. [Tq| ), 

with slopes = dXi/dXp^^^^ a'^ = daQ/dXp^^^^ a'^ = dai/dXp^^^. In this formalism we have only four 
adjustable parameters Xc, Aq, and we can obtain directly an uncertainty on according to a 

given hadronic model. In the following computations, we use as derivative for the fourth parameter: 

dN dN , dN , dN 

^ 3Ai dao dai 

For a set of measurements the logarithm of the likelihood may be written, omitting the constant terms, 
as: 

^(p) = k* V’i; p)) - N{ri, 'ipi; p)) (16) 

i 

Maximizing the likelihood gives an estimator of p with an error (covariance) matrix C = W~^ where the 
weight matrix W is defined as: 


Wjk = - 


dpjdpk 


Uj dNj dNj 

, \ V? dpj dpk 


(1-^ 


d^Ni 


Ni dpjdpk 


(17) 


where the index i for N and its derivatives means “at position In average rii = so we obtain a 

simple average expression for the elements of W: 


^ ^ 1 dNi dNi d{\nNi) d{\nNi 

^ Ni dpj dpk ^ 


dpi 


dpk 


(18) 


Various simulations have shown that the weight matrix does not vary strongly in different realizations, so 
we take the inverse of IF as a good estimator of the error matrix we can obtain in real measurements 


A direct consequence of Eq 18 is the scaling property of the errors: the derivatives d{\nNi)/dpj depend 


only on the position, so if V(r, %p) is multiplied by a global factor F, the weight matrix is multiplied by F, 
and the errors are divided by ^/F. As a consequence, for a given geometrical configuration (^, Bt and array 
of detectors), the errors are approximately proportional to l/y^Fprim- A scaling with the spacing I may be 
obtained in the approximation of a dense array, if the summation in Eqj^ may be replaced by an integral: 
this is reasonable here because ln(V(r, is a smooth function over the front plane, except at the origin. As 
a result, IF is proportional to l/(^^ cos6>), so the errors are proportional to iVcosO^ that is to the inverse of 
the square root of the density of detectors in projection onto the front plane. 


6.3. Precision on the indirect measurement of the depth of maximum 

Using the values Ai, A'^, ao, o^q, <ai, ol'i obtained from our sample of simulations of Proton and iron showers 
at 1 EeV at different values of 0 and Bt, we can evaluate using Eqj^the errors we can expect on an indirect 


measurement of the results for the standard array with r^ef = 1000 m are plotted on Fig. 14 for 


Ao = 0, that is V(rref,7r/4) = 1 (normalized size of the muon component). The errors on ao, ai and Ai 
depend little on Bt and moderately on 6>, while the error on from the compound model decreases 

strongly, as expected, with increasing Bt. The dependence on 0 is more complicated: for Bt > 30/iT, the 
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precision on is dominated by the dependence of o^o, so the error decreases with increasing 0] at low 

0 and the dependence of Ai is the most important one and may favour a low value of 0. If the errors 
are computed for a fixed energy, the dependence on 0 is compensated by the decrease of N with 0; using 
the scaling properties, the result may be summarized in the following way: for the highest values of the 
error on (in g/cm^) is of the order of 20W\/EA g/cm^, with I in km, E in EeV, A in m^. As a result, 

when using a given hadronic model, an event by event identification may be envisaged if the error on X!^^^ 
is about 50 g/cw? or less, that is, for the standard array, at energies of the order of 5 EeV or more, if the 
local geomagnetic field is at least 40 /iT. 
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Figure 14: Errors on the shape parameters ao, oti and Ai and induced on X^ax as functions of 9 and Bt. Top: with a 
normalized size of the muon component (1 hit in average in a detector at 1000 m from core, at 45 degrees from the magnetic 
field in the transverse plane); bottom: scaled at 1 EeV, with an effective area of 10 m^ per detector). 


If the muon detector array is completed by an independent measurement of 
tightly correlated), comparing Eig. 
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(or Xmax, which is 
and Eigl4 indicates in which conditions the shape parameters can 


provide some discrimination between hadronic models. 


7. Application to possible detectors 

The results found here may be used at different levels of demanding observations, that is also at different 
levels of dependence on models. 

• Assuming a reliable model of hadronic interaction: an array of muon detectors can obtain an average 
value of X^^^ at lower energies and possibly an event-by-event determination at high energy. As an 
internal check, the model has to reproduce the observed dependence on 0. 

• Using an external calibration, that is taking the average X^ax from another experiment in similar 
conditions of zenith angle, and correcting for a different atmosphere profile if needed. The hadronic 
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model has to give the right average value of if the model is acceptable, the spectrum of 

may be exploited, especially on an event-by-event basis at highest energies. This option may be uneasy, 
because it requires a combination of different experiments. 

• Building a hybrid detector including a muon array and a longitudinal profile detector covering a few 
ten kilometers in at least one direction from the muon array, a self-calibrated measurement could be 
achieved. The profile detector could be made of single fluorescence eyes with a wide field of view as 
proposed in [20]: for nearly horizontal showers passing above such eyes, the time profile of the received 
light gives a good measurement of the portion of profile within the field of view. A possible layout of 
a hybrid detector is proposed in Fig. Any other detector able to see the profile from the side could 
replace the fluorescence detection with the advantage of a larger duty cycle. 



Figure 15: A possible layout for a hybrid detector. Left part: muon array: vertical detectors with two orientations to obtain a 
roughly isotropic sensitivity to horizontal showers; the red ellipses are the contour on ground at 1000 m from the shower axes, 
for showers at ^ = 70 deg. Right part: fluorescence eyes or other longitudinal profile detectors; the blue curve is a typical 
longituflnal profile for such showers. 


The intensity and the inclination of the magnetic field depends on the location on the surface on the 
Earth. It has a maximum of about 65 jaT in the North of Asia or in the antarctic region South of Australia, 
where its direction is nearly vertical, providing a large value of Bf for nearly horizontal incidence, whatever 
the azimuth angle. 

Another important parameter is the altitude. Higher altitude is preferable for higher zenith angle, 
because the attenuation of the muon component may prevent an efficient detection after a long path; in 
the case of a hybrid detector, another advantage is also to reduce the distance between the ground and the 
maximum of the longitudinal profile. 

In this study, we have assumed an horizontal ground, but this is not an experimental requirement. The 
slope of the ground may increase the aperture for a restricted azimuthal region, and an elongated shape of 
the array, with a larger spacing of the detectors in that direction, may be the best option. 

In any configuration, the profile of the atmosphere needs to be known precisely, because it may be the 
main source of systematic errors. In case of large diurnal and/or seasonal variations, it should be carefully 
monitored. 

8. Summary 

Using a semi-analytic calculation of the density of muons in inclined atmospheric showers, in the approx¬ 
imation of small angular deviations, we have shown that the geomagnetic deflection provides an indirect 
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measurement of the depth of maximal production of muons, which is an indicator of the nature of 

the primary particle, tightly correlated to the usual parameter Xmax- This measurement may be performed 
using an universal parametrization of the muon density as a function of the distance to the shower axis and 
the azimuthal angle in projection onto the front plane (perpendicular to the shower axis at ground level). 
A significant dependence on was found for three parameters Ai, ai and a 2 ^ especially the second one 

which describes the amplitude of the quadrupolar distortion of the muon density in the front plane due to 
the transverse component Bf of the magnetic field. The calculation of the density was successfully compared 
to detailed simulations of the muon propagation in some showers. 

The parametrization may be inaccurate for very strong distortions, that is for highest values of 0 and/or 
Bf. in such cases the sensitivity to X^ax is certainly strong and a more accurate simulation of trajectories 
and a more complex parametrization may be needed to fully exploit the observations. However at large 0 
the muonic flux is attenuated and more difficult to measure at ground level; the potential gain on primary 
identification using events with 6> > 80 needs a specific study (maybe with a more elaborated parametriza¬ 
tion), beyond the scope of this paper. 

A precise measurement on the direction is needed as the dependence of the parameters on 0 is important, 
and increases with 6>; fortunately the front of inclined showers is thin, and the thickness decreases with 6>, 
so a precise measurement of the arrival time of the muons (at the level of 10 ns or better) can determine 0 
with enough precision to extract the dependence of the parameters on 

In practice the density is measured on the ground instead of the front plane, so the longitudinal evolution 
of the shower and the divergence of muons from the axis result in a forward/backward asymmetry of the 
observed density: this effect is essentially a dipolar distortion which may be evaluated and corrected for; in 
any case it may be disentangled from the magnetic distortion. The linear relation between the distortion 
parameters and X^^^^ is slightly shifted by a quantity depending on the orientation of the transverse field 
within the front plane, but the slope is not modified. A correction for asymmetries of the front at ground 
level may be needed to suppress a possible bias on 0] detailed simulations can estimate such a bias. 

The variations of the atmospheric density profile at low altitude (below the region of production of the 
muons) impact mainly the energy loss, so the usual variations of a few percent do not modify too much the 
density at ground, and they are easy to monitor. Variations in upper atmosphere are in principle not large, 
but they may affect the development of the hadronic cascade and modify the altitude of production of the 
muons and their energy spectrum: this may change both the position of X^^^^ with respect to Xmax, and 
the relation between the distortion parameters and which depends on the distance from production 

to ground. So for a given site, the profile of the upper atmosphere and its possible variations should be known. 

We have proposed some possible scenarios of application of the mass sensitive parameters measured 
with the method developped in this paper. Especially, a hybrid detector able to measure simultaneously 
the longitudinal profile and the magnetic distortion of muons in horizontal showers, can provide strong 
constraints on the models of hadronic interactions. This option could be valuable to step forward in both 
the hadronic modelling and in the determination of the composition of UHECRs. 
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Appendix: analytic computation of the ground density with an exponential atmosphere 

To simplify the equations we use the coordinates of the shower frame (z axis along the shower axis), and 
we omit the si subscripts. The depth is described by the function: 

X(z) = Xa exp(-z/L) (19) 


where Xa and L account for the slant of the axis. 

The muons which reach the ground are ultrarelativistic, except possibly in the very end of the trajectory. 
So we use E = p everywher^ Their direction remains close to the z axis so we can use small angle 
approximations. 

Let us consider a muon of mass m and lifetime r, injected at Zi with energy Ei. When going down from 
z to z — d^:, we have to account for: 

• the decay probability: {—dz)/{ ctE/ m) = {—(lz)/{\E) with X^crjm. 

• the energy loss: dE = —edX = —eXajE exp{—z/L)dz. Introducing the energy extrapolated 
backwards to infinity Eoo = Ei d- sXajE exp(— 2 :^/!/), we can write: E{z) = Ei — eAX = E^o — 
eXaexp{-z/L) 

• the magnetic deflection (in the direction perpendicular to the transverse field 5^): duj = P/E{z)dz 
with P = cBt if E is expressed in eV. This elementary deflection produces a deviation zduj when 
extrapolated to z = 0. 

• multiple scattering: a random angular deflection with a variance dr]^ = {Ems/E{z))^ dX/X^ad = 
— {EmslE{z))‘^e~^^^dzl{LXrad) in both transverse directions, where X^ad is the radiation length of 
air and Ej^s = 13.6 MeV. 

The cumulative effects at ground level {z = 0, energy Eg^ = E^^ — eXa) are then: 


• survival probability P{z)\ 


1 dP 

~P~^ 

InP(O) 


\{E^-eXae--/i^) 


In 


XEoo V ^ 00 
L , fXiE„ 


\E^ 


In 


Jgj> 


x„, Ei 


d(lnP) 


L d(e^/-^) 

XE^ - eXa/Eoc, 


XEq^ \ ^ 9 ^ 


• angular deviation from Zi to 0: 


dz 


UJ = 


E^-XaSe-^/L- 

pL ^( e-^l^-eXalE, 
Eoo V 1~£W/-Eoo 
Ei 


pL 




Eoo lo - eXajEoo 

pL ( (l-e-^^/^eXa/Eo 

= -ET- Zi/L + In 


Eo 


1 — ^Xa/Eo 


E^ 


Lin- 


E, 


gr 


^For convenience we omit c when using the momentum p and the mass m, that is, we write p for pc and m for mc^; however, 
we keep cr in the expression of the decay length. 
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• position deviation from Zi to 0; 


5 = p 


f 


dz 


Eoo - XaSe-^l^ E, 


= ^E,{a,Zi/L) 


with 


a = —— and hi 

Eqo 




udu 


— ae 


• variance of deviation by multiple scattering: we make a summation over independent angular deviations 
in atmosphere slices: 


ang 


Ei 




Xae-^/^dz 

LE{zfXrad 
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\ -^CJO / Xj-Q^d Jq 
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Ei 


Jgj. 


Ei 


variance of the deviation: we make a summation over the contributions of the independent angular 
deviations to the position at z = 


= Ei 


Xge-^/^z^dz ^ fEr^V Xg r 
Jo LE{z)‘^Xrad \ Eoo ) XradL Jq 


Zi g-z/L^2 


XradL Jo (l-ae --/^)2 


E, 


ms \ XgL . Zi. 

V-^2(a,y) 

-^oo / -^rad ^ 


with 


F 2 {a,C)= [ 

Jo 


^ e du 

(1 — ae~^Y 


The function E 2 may be expressed with Ei and Ei may be expressed trough the special function Li 2 
(dilogarithm): 


2 Ei{aX) 


e 


^2(0^, C) = 

a a[l — ae 

(-2 

Ei{a, C) = ^ + Cln(l - ae~^) + Li2((a) - lj\2{ae~^) 

ln(l — u) 


2 

k 


ri.W = gF'-i 


du 


The functions Ei and E 2 are plotted in Fig. 16 for different values of a between 0 and 1 (if a > 1 the muon 
stops before reaching the ground). 
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Figure 16: Functions Fi{x)lx^ (left) and F 2 {x) (right) used to express the magnetic deviation and the multiple scattering 
dispersion at ground level. From bottom to top: a = 0.1, 0.2, • • •, 1 
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